1 Setup

Set the working dir

setwd("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/Conditions")

Load libraries

suppressPackageStartupMessages({
    library(data.table)
    library(DESeq2)
    library(gplots)
    library(here)
    library(hyperSpec)
    library(limma)
    library(LSD)
    library(magrittr)
    library(matrixStats)
    library(parallel)
    library(pander)
    library(plotly)
    library(plyr)
    library(RColorBrewer)
    library(scatterplot3d)
    library(tidyverse)
    library(tximport)
    library(VennDiagram)
    library(vsn)
})

Source some helper functions

suppressPackageStartupMessages({
    source("~/Git/UPSCb/UPSCb-common/src/R/featureSelection.R")
    source("~/Git/UPSCb/UPSCb-common/src/R/gopher.R")
    source("~/Git/UPSCb/UPSCb-common/src/R/plot.multidensity.R")
    source("~/Git/UPSCb/UPSCb-common/src/R/volcanoPlot.R")
})

Cutoff Different from Schurch et al., RNA, 2016 (LFC at 1 instead of 0,5)

lfc <- 1
FDR <- 0.01

Create palettes

pal <- c(brewer.pal(8,"Dark2"),1)
pal2 <- brewer.pal(9,"Paired") #require package RColorBrewer
cols <- rainbow(17)
hpal <- colorRampPalette(c("blue","white","red"))(100)

mar <- par("mar")

2 Expression data

2.1 Loading of Fold changes and adjusted p-values

Read the expression data for 6 hours

T6_NP_DMSO <- read.csv("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/Conditions/Conditions_T6_NP_DMSO_vs_T6_NPS_DMSO.csv",sep=";",dec=",")
T6_NS_DMSO <- read.csv("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/Conditions/Conditions_T6_NS_DMSO_vs_T6_NPS_DMSO.csv",sep=";",dec=",")
T6_NPS_AZD <- read.csv("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/Conditions/Conditions_T6_NPS_AZD_vs_T6_NPS_DMSO.csv",sep=";",dec=",")
T6_NP_AZD <- read.csv("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/Conditions/Conditions_T6_NP_AZD_vs_T6_NPS_DMSO.csv",sep=";",dec=",")
T6_NS_AZD <- read.csv("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/Conditions/Conditions_T6_NS_AZD_vs_T6_NPS_DMSO.csv",sep=";",dec=",")

Read the expression data for 24 hours

T24_NP_DMSO <- read.csv("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/Conditions/Conditions_T24_NP_DMSO_vs_T24_NPS_DMSO.csv",sep=";",dec=",")
T24_NS_DMSO <- read.csv("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/Conditions/Conditions_T24_NS_DMSO_vs_T24_NPS_DMSO.csv",sep=";",dec=",")
T24_NPS_AZD <- read.csv("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/Conditions/Conditions_T24_NPS_AZD_vs_T24_NPS_DMSO.csv",sep=";",dec=",")
T24_NP_AZD <- read.csv("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/Conditions/Conditions_T24_NP_AZD_vs_T24_NPS_DMSO.csv",sep=";",dec=",")
T24_NS_AZD <- read.csv("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/Conditions/Conditions_T24_NS_AZD_vs_T24_NPS_DMSO.csv",sep=";",dec=",")

Read Wouter’s data

Seedling <- read.csv("/mnt/picea/projects/arabidopsis/jhanson/arabidopsis-nutrition-TOR/Wouter_short.csv",sep=";",dec=",")
Seedling <- cbind(Seedling["x"],Seedling["Coef.acol...col"],Seedling["p.value.adj.acol...col"])
colnames(Seedling) <- c("rownames.res.","res.log2FoldChange","res.padj")

2.2 Apply filters for the DEGs

2.2.1 For the 6hrs timepoint

res <- T6_NP_DMSO
cutoff1 <- res$res.log2FoldChange >= lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
cutoff2 <- res$res.log2FoldChange <= -lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
SucLow6hrs <- res$rownames.res.[cutoff2]
SucHigh6hrs <- res$rownames.res.[cutoff1]

res <- T6_NS_DMSO
cutoff1 <- res$res.log2FoldChange >= lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
cutoff2 <- res$res.log2FoldChange <= -lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
PiLow6hrs <- res$rownames.res.[cutoff2]
PiHigh6hrs <- res$rownames.res.[cutoff1]

res <- T6_NPS_AZD
cutoff1 <- res$res.log2FoldChange >= lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
cutoff2 <- res$res.log2FoldChange <= -lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
AzdLow6hrs <- res$rownames.res.[cutoff2]
AzdHigh6hrs <- res$rownames.res.[cutoff1]

res <- T6_NP_AZD
cutoff1 <- res$res.log2FoldChange >= lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
cutoff2 <- res$res.log2FoldChange <= -lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
AzdNPLow6hrs <- res$rownames.res.[cutoff2]
AzdNPHigh6hrs <- res$rownames.res.[cutoff1]

res <- T6_NS_AZD
cutoff1 <- res$res.log2FoldChange >= lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
cutoff2 <- res$res.log2FoldChange <= -lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
AzdNSLow6hrs <- res$rownames.res.[cutoff2]
AzdNSHigh6hrs <- res$rownames.res.[cutoff1]

2.2.2 For the 24hrs timepoint

res <- T24_NP_DMSO
cutoff1 <- res$res.log2FoldChange >= lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
cutoff2 <- res$res.log2FoldChange <= -lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
SucLow24hrs <- res$rownames.res.[cutoff2]
SucHigh24hrs <- res$rownames.res.[cutoff1]

res <- T24_NS_DMSO
cutoff1 <- res$res.log2FoldChange >= lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
cutoff2 <- res$res.log2FoldChange <= -lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
PiLow24hrs <- res$rownames.res.[cutoff2]
PiHigh24hrs <- res$rownames.res.[cutoff1]

res <- T24_NPS_AZD
cutoff1 <- res$res.log2FoldChange >= lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
cutoff2 <- res$res.log2FoldChange <= -lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
AzdLow24hrs <- res$rownames.res.[cutoff2]
AzdHigh24hrs <- res$rownames.res.[cutoff1]

res <- T24_NP_AZD
cutoff1 <- res$res.log2FoldChange >= lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
cutoff2 <- res$res.log2FoldChange <= -lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
AzdNPLow24hrs <- res$rownames.res.[cutoff2]
AzdNPHigh24hrs <- res$rownames.res.[cutoff1]

res <- T24_NS_AZD
cutoff1 <- res$res.log2FoldChange >= lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
cutoff2 <- res$res.log2FoldChange <= -lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
AzdNSLow24hrs <- res$rownames.res.[cutoff2]
AzdNSHigh24hrs <- res$rownames.res.[cutoff1]

2.2.3 For Wouter’s data

res <- as.data.frame(Seedling)
cutoff1 <- res$res.log2FoldChange >= lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
cutoff2 <- res$res.log2FoldChange <= -lfc & ! is.na(res$res.padj) & res$res.padj <= FDR
Seedling_Low <- res$rownames.res.[cutoff2]
Seedling_High <- res$rownames.res.[cutoff1]

3 GO analysis

3.1 For the 6hr timepoint

GO <- gopher(SucLow6hrs,url="athaliana", alpha=2)

## Loading required package: jsonlite

## 
## Attaching package: 'jsonlite'

## The following object is masked from 'package:purrr':
## 
##     flatten

## No enrichments found in task: mapman

GO_SucLow6hrs <- GO$go

GO <- gopher(SucHigh6hrs,url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_SucHigh6hrs <- GO$go

GO <- gopher(PiLow6hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_PiLow6hrs <- GO$go

GO <- gopher(PiHigh6hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_PiHigh6hrs <- GO$go

GO <- gopher(AzdLow6hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_AzdLow6hrs <- GO$go

GO <- gopher(AzdHigh6hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_AzdHigh6hrs <- GO$go

GO <- gopher(AzdNPLow6hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_AzdNPLow6hrs <- GO$go

GO <- gopher(AzdNPHigh6hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_AzdNPHigh6hrs <- GO$go

GO <- gopher(AzdNSLow6hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_AzdNSLow6hrs <- GO$go

GO <- gopher(AzdNSHigh6hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_AzdNSHigh6hrs <- GO$go

3.2 For the 24hr timepoint

GO <- gopher(SucLow24hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_SucLow24hrs <- GO$go

GO <- gopher(SucHigh24hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_SucHigh24hrs <- GO$go

GO <- gopher(PiLow24hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_PiLow24hrs <- GO$go

GO <- gopher(PiHigh24hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_PiHigh24hrs <- GO$go

GO <- gopher(AzdLow24hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_AzdLow24hrs <- GO$go

GO <- gopher(AzdHigh24hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_AzdHigh24hrs <- GO$go

GO <- gopher(AzdNPLow24hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_AzdNPLow24hrs <- GO$go

GO <- gopher(AzdNPHigh24hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_AzdNPHigh24hrs <- GO$go

GO <- gopher(AzdNSLow24hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_AzdNSLow24hrs <- GO$go

GO <- gopher(AzdNSHigh24hrs, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_AzdNSHigh24hrs <- GO$go

3.3 For Wouter’s data

GO <- gopher(Seedling_Low, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_Seedling_Low <- GO$go

GO <- gopher(Seedling_High, url="athaliana", alpha=2)

## No enrichments found in task: mapman

GO_Seedling_High <- GO$go

3.4 Add the experiment name in the GO tables

GO_AzdLow6hrs <- cbind(Description='GO_Azd6hrs', GO_AzdLow6hrs)
GO_AzdHigh6hrs <- cbind(Description='GO_Azd6hrs', GO_AzdHigh6hrs)
GO_SucHigh6hrs <- cbind(Description='GO_Suc6hrs', GO_SucHigh6hrs)
GO_SucLow6hrs <- cbind(Description='GO_Suc6hrs', GO_SucLow6hrs)
GO_PiHigh6hrs <- cbind(Description='GO_Pi6hrs', GO_PiHigh6hrs)
GO_PiLow6hrs <- cbind(Description='GO_Pi6hrs', GO_PiLow6hrs)
GO_AzdNPHigh6hrs <- cbind(Description='GO_AzdNP6hrs', GO_AzdNPHigh6hrs)
GO_AzdNPLow6hrs <- cbind(Description='GO_AzdNP6hrs', GO_AzdNPLow6hrs)
GO_AzdNSHigh6hrs <- cbind(Description='GO_AzdNS6hrs', GO_AzdNSHigh6hrs)
GO_AzdNSLow6hrs <- cbind(Description='GO_AzdNS6hrs', GO_AzdNSLow6hrs)
GO_AzdHigh24hrs <- cbind(Description='GO_Azd24hrs', GO_AzdHigh24hrs)
GO_AzdLow24hrs <- cbind(Description='GO_Azd24hrs', GO_AzdLow24hrs)
GO_SucHigh24hrs <- cbind(Description='GO_Suc24hrs', GO_SucHigh24hrs)
GO_SucLow24hrs <- cbind(Description='GO_Suc24hrs', GO_SucLow24hrs)
GO_PiHigh24hrs <- cbind(Description='GO_Pi24hrs', GO_PiHigh24hrs)
GO_PiLow24hrs <- cbind(Description='GO_Pi24hrs', GO_PiLow24hrs)
GO_AzdNPHigh24hrs <- cbind(Description='GO_AzdNP24hrs', GO_AzdNPHigh24hrs)
GO_AzdNPLow24hrs <- cbind(Description='GO_AzdNP24hrs', GO_AzdNPLow24hrs)
GO_AzdNSHigh24hrs <- cbind(Description='GO_AzdNS24hrs', GO_AzdNSHigh24hrs)
GO_AzdNSLow24hrs <- cbind(Description='GO_AzdNS24hrs', GO_AzdNSLow24hrs)
GO_Seedling_High <- cbind(Description='GO_Seedling', GO_Seedling_High)
GO_Seedling_Low <- cbind(Description='GO_Seedling', GO_Seedling_Low)

4 Combine the GO tables

4.1 Prepare the tables

combined_up <- rbind.fill(GO_AzdHigh6hrs, GO_AzdNPHigh6hrs, GO_AzdNSHigh6hrs, GO_PiHigh6hrs, GO_SucHigh6hrs,GO_AzdHigh24hrs, GO_AzdNPHigh24hrs, GO_AzdNSHigh24hrs, GO_PiHigh24hrs, GO_SucHigh24hrs,GO_Seedling_High)
combined_down <- rbind.fill(GO_AzdLow6hrs,GO_SucLow6hrs,GO_PiLow6hrs,GO_AzdNPLow6hrs,GO_AzdNSLow6hrs,GO_AzdLow24hrs,GO_SucLow24hrs,GO_PiLow24hrs,GO_AzdNPLow24hrs,GO_AzdNSLow24hrs,GO_Seedling_Low)

4.2 Keep only the biological processes

combined_up <- combined_up[combined_up$namespace == "BP",]
combined_down <- combined_down[combined_down$namespace == "BP",]

4.3 Combining Up and Down-regulated genes together

combined_up[,"Direction"]  <- "Induced"
combined_down[,"Direction"]  <- "Repressed"

combined_all_up_and_down <- rbind(combined_up, combined_down)

Export the obtained data.frame containing also the non-significant GOs

write.table(combined_all_up_and_down,file="GOs_also_NON_sgnificant_GOs.txt", quote = F ,sep = "\t", row.names=FALSE)

4.4 Removing the root terms

exclude <- c("GO:0008150", "GO:0008283", "GO:0032502", "GO:0044848", "GO:0032501", "GO:0098743", "GO:0008152", "GO:0051704", "GO:0007610", "GO:0000003", "GO:0009987",
             'GO:0006791', 'GO:0001906', 'GO:0002376', 'GO:0006794', 'GO:0009758', 'GO:0015976', 'GO:0019740', 'GO:0022414', 
             'GO:0022610', 'GO:0023052', 'GO:0040007', 'GO:0040011', 'GO:0043473', 'GO:0048511', 'GO:0048518', 'GO:0048519', 'GO:0050789', 'GO:0050896', 'GO:0051179', 'GO:0065007', 'GO:0071840', 'GO:0098754', 'GO:0110148')
combined_all_up_and_down <- combined_all_up_and_down[! combined_all_up_and_down$id %in% exclude,]

# ' ## Change the p-values to change all the ones that are non-significant
combined_all_up_and_down$padj[combined_all_up_and_down$padj > 0.05] <- NA #Change this value as wished (normally 0.05)

5 Create a list of interesting GO

Create a list of cell cycle related GO enriched in Repressed genes after 24hrs of Pi starvation

C_Cycle <- c("GO:0007049","GO:0051726","GO:0044770","GO:0007051","GO:0007059","GO:0042549","GO:0051301","GO:0051304",
             "GO:0009061","GO:0007017","GO:0046490","GO:0006720","GO:0032196")

Pyrimidine-containing compound metabolism

Pyrimidine <- c("GO:0072527",'GO:0006304','GO:0055086','GO:0018193','GO:0016072','GO:0090501','GO:0006354',
                'GO:1901564','GO:1901566','GO:0006139','GO:0000375','GO:0043603','GO:0032446','GO:0019318','GO:0006259',
                'GO:0006260','GO:0090305','GO:0006518','GO:0006553','GO:0006090','GO:0006089',
                'GO:0008213','GO:0006396','GO:0034641','GO:0005991','GO:0009451','GO:0034660')

DNA replication

Replication <- c('GO:0006260','GO:0006304','GO:0090066','GO:0055086','GO:0018193','GO:0048519','GO:0006793',
                 'GO:0031323','GO:0045814','GO:0006270','GO:0006275','GO:0018205','GO:0016072','GO:0044260',
                 'GO:0090501','GO:0009069','GO:0050789','GO:1901576','GO:0046483','GO:1901564','GO:0051338',
                 'GO:0006139','GO:0043603','GO:0032446','GO:0006725','GO:0019318','GO:0006259',
                 'GO:0010498','GO:0006518','GO:0010467','GO:0033559','GO:0006090','GO:0043412','GO:0008213',
                 'GO:0060249','GO:0040029','GO:0009892','GO:0040008','GO:0072527','GO:0034641',
                 'GO:0005991','GO:0009451','GO:0034660')

Lipopolysaccharide metabolism (repressed by sugar starvation after 24 hrs)

LPS <- c("GO:1903509","GO:1901137","GO:0006643","GO:0006629")

Response to heat (repressed by sugar starvation after 24 hrs)

heat <- c("GO:0009408","GO:0071496","GO:0009642","GO:0016036","GO:0006950","GO:0006979","GO:1901700","GO:0002215")

ncRNA (induced by 24hrs of sugar starvation)

ncRNA <- c('GO:0034660','GO:0042537','GO:0009768','GO:0055070','GO:0043900','GO:0051246','GO:0018198','GO:0017014',
           'GO:0016054','GO:0016070','GO:0016072','GO:0042592','GO:0022900','GO:0055114','GO:0006081','GO:0006082',
           'GO:0008219','GO:0009069','GO:0009308','GO:0001505','GO:0006396','GO:1901605','GO:0030258','GO:0044283',
           'GO:0009850','GO:0044281','GO:0005991','GO:0005997','GO:0042430','GO:0051174','GO:0065008','GO:1901657','GO:0006766')


myList <- list(C_Cycle,Pyrimidine,Replication,LPS, heat, ncRNA)

5.1 List of GO groups overrepresented among genes induced after 24hrs of Suc starvation

aging <- c("GO:0007568","GO:0090558","GO:0048366","GO:0010618","GO:0007275")
response_to_red_light <- c("GO:0010114","GO:1901700","GO:0009743","GO:0006290","GO:0080021","GO:0017085","GO:0043617","GO:0034285","GO:0009637","GO:0009646","GO:0009611","GO:0010033")
photosynthesis <- "GO:0015979"
response_to_stimulus <- "GO:0050896"
localization <- "GO:0051179"
multi_organism_process <- "GO:0051704"
biological_regulation <- "GO:0065007"
catabolism <- "GO:0009056"
plastid_organization <- c("GO:0009657","GO:0043933","GO:0072657","GO:0022613","GO:0008643","GO:0044085","GO:0010256")
isoprenoid_transport <- c("GO:0046864","GO:0055085","GO:0006811","GO:0015850","GO:0032879","GO:0015853","GO:0072511","GO:0010232","GO:0015672","GO:0015718")
drug_metabolism <- "GO:0017144"
multidimensional_cell_growth <- "GO:0009825"
generation_of_precursor_metabolites_and_energy <- "GO:0006091"
secondary_metabolism <- c("GO:0019748","GO:0015977")
ROS_metabolism <- "GO:0072593"
sulfur_compound_biosynthesis <- "GO:0044272"
dephosphorylation <- c("GO:0016311","GO:0006733","GO:0006739")
sulfur_compound_metabolism <- "GO:0006790"
ncRNA_metabolism <- c("GO:0034660","GO:0042537","GO:0009768","GO:0055070","GO:0043900","GO:0006520","GO:0051246","GO:0018198","GO:0017014","GO:0048583","GO:0016070","GO:0016072","GO:0042592","GO:0033559","GO:0022900","GO:0055114","GO:0006081","GO:0006082","GO:0008219","GO:0009311","GO:0009069","GO:0009308","GO:0001505","GO:0006396","GO:0009850","GO:0044283","GO:0044281","GO:0009607","GO:0005991","GO:0044093","GO:0005997","GO:0042430","GO:0051174","GO:0065008","GO:0042221","GO:0006766")
cofactor_metabolism <- c("GO:0051186","GO:0006793")
carbohydrate_metabolism <- "GO:0005975"
ureide_metabolism <- "GO:0010135"
oxidative_photosynthetic_carbon_pathway <- "GO:0009854"


myList <- list(
    aging = aging,
    response_to_red_light= response_to_red_light,
    photosynthesis=photosynthesis,
    catabolism=catabolism,
    plastid_organization=plastid_organization,
    isoprenoid_transport=isoprenoid_transport,
    drug_metabolism=drug_metabolism,
    multidimensional_cell_growth=multidimensional_cell_growth,
    generation_of_precursor_metabolites_and_energy=generation_of_precursor_metabolites_and_energy,
    secondary_metabolism=secondary_metabolism,
    ROS_metabolism=ROS_metabolism,
    sulfur_compound_biosynthesis=sulfur_compound_biosynthesis,
    dephosphorylation=dephosphorylation,
    sulfur_compound_metabolism=sulfur_compound_metabolism,
    ncRNA_metabolism=ncRNA_metabolism,
    cofactor_metabolism=cofactor_metabolism,
    carbohydrate_metabolism=carbohydrate_metabolism,
    ureide_metabolism=ureide_metabolism,
    oxidative_photosynthetic_carbon_pathway=oxidative_photosynthetic_carbon_pathway)

5.2 List of GO groups overrepresented under the different conditions

LowSuc6hr <- GO_SucLow6hrs$id[GO_SucLow6hrs$padj <= 0.05]
LowPi6hr <- GO_PiLow6hrs$id[GO_PiLow6hrs$padj <= 0.05]
LowAzd6hr <- GO_AzdLow6hrs$id[GO_AzdLow6hrs$padj <= 0.05]
LowAzdSuc6hr <- GO_AzdNPLow6hrs$id[GO_AzdNPLow6hrs$padj <= 0.05]
LowAzdPi6hr <- GO_AzdNSLow6hrs$id[GO_AzdNSLow6hrs$padj <= 0.05]

HighSuc6hr <- GO_SucHigh6hrs$id[GO_SucHigh6hrs$padj <= 0.05]
HighPi6hr <- GO_PiHigh6hrs$id[GO_PiHigh6hrs$padj <= 0.05]
HighAzd6hr <- GO_AzdHigh6hrs$id[GO_AzdHigh6hrs$padj <= 0.05]
HighAzdSuc6hr <- GO_AzdNPHigh6hrs$id[GO_AzdNPHigh6hrs$padj <= 0.05]
HighAzdPi6hr <- GO_AzdNSHigh6hrs$id[GO_AzdNSHigh6hrs$padj <= 0.05]

LowSuc24hr <- GO_SucLow24hrs$id[GO_SucLow24hrs$padj <= 0.05]
LowPi24hr <- GO_PiLow24hrs$id[GO_PiLow24hrs$padj <= 0.05]
LowAzd24hr <- GO_AzdLow24hrs$id[GO_AzdLow24hrs$padj <= 0.05]
LowAzdSuc24hr <- GO_AzdNPLow24hrs$id[GO_AzdNPLow24hrs$padj <= 0.05]
LowAzdPi24hr <- GO_AzdNSLow24hrs$id[GO_AzdNSLow24hrs$padj <= 0.05]

HighSuc24hr <- GO_SucHigh24hrs$id[GO_SucHigh24hrs$padj <= 0.05]
HighPi24hr <- GO_PiHigh24hrs$id[GO_PiHigh24hrs$padj <= 0.05]
HighAzd24hr <- GO_AzdHigh24hrs$id[GO_AzdHigh24hrs$padj <= 0.05]
HighAzdSuc24hr <- GO_AzdNPHigh24hrs$id[GO_AzdNPHigh24hrs$padj <= 0.05]
HighAzdPi24hr <- GO_AzdNSHigh24hrs$id[GO_AzdNSHigh24hrs$padj <= 0.05]

LowSeedling <- GO_Seedling_Low$id[GO_Seedling_Low$padj <= 0.05]
HighSeedling <- GO_Seedling_High$id[GO_Seedling_High$padj <= 0.05]

myList <- list(LowSuc6hr = LowSuc6hr,
               LowPi6hr = LowPi6hr,
               LowAzd6hr= LowAzd6hr,
               LowAzdSuc6hr= LowAzdSuc6hr,
               LowAzdPi6hr= LowAzdPi6hr,
               HighSuc6hr= HighSuc6hr,
               HighPi6hr= HighPi6hr,
               HighAzd6hr= HighAzd6hr,
               HighAzdSuc6hr= HighAzdSuc6hr,
               HighAzdPi6hr= HighAzdPi6hr,
               LowSuc24hr= LowSuc24hr,
               LowPi24hr= LowPi24hr,
               LowAzd24hr= LowAzd24hr,
               LowAzdSuc24hr= LowAzdSuc24hr,
               LowAzdPi24hr= LowAzdPi24hr,
               HighSuc24hr= HighSuc24hr,
               HighPi24hr= HighPi24hr,
               HighAzd24hr= HighAzd24hr,
               HighAzdSuc24hr= HighAzdSuc24hr,
               HighAzdPi24hr= HighAzdPi24hr,
               LowSeedling = LowSeedling,
               HighSeedling=HighSeedling
)

6 Create a sublist of GO with only the GOs of interest

go <- combined_all_up_and_down[combined_all_up_and_down$id %in% myList,]
go <- lapply(myList,function(x)
{
    combined_all_up_and_down[combined_all_up_and_down$id %in% x,]
})

7 Make the bubble plot

lapply(names(go),function(nam,l)
{
    x <- l[[nam]]
    ggplot(x,aes(x=Description,y=id))+ geom_point(aes(size=npat, color=-log10(padj)))+ 
        facet_wrap(~Direction, nrow=1) +  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
        ggtitle(nam)+scale_shape_discrete(solid=F)+
        scale_colour_gradientn(colours=c("red2", "purple2", "blue2")) +
        theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 1)) 
}, go)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

8 Session Info

## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] jsonlite_1.6                vsn_3.54.0                 
##  [3] VennDiagram_1.6.20          futile.logger_1.4.3        
##  [5] tximport_1.14.0             forcats_0.4.0              
##  [7] stringr_1.4.0               dplyr_0.8.3                
##  [9] purrr_0.3.3                 readr_1.3.1                
## [11] tidyr_1.0.0                 tibble_2.1.3               
## [13] tidyverse_1.3.0             scatterplot3d_0.3-41       
## [15] RColorBrewer_1.1-2          plyr_1.8.4                 
## [17] plotly_4.9.1                pander_0.6.3               
## [19] magrittr_1.5                LSD_4.0-0                  
## [21] limma_3.42.0                hyperSpec_0.99-20180627    
## [23] ggplot2_3.2.1               lattice_0.20-38            
## [25] here_0.1                    gplots_3.0.1.1             
## [27] DESeq2_1.26.0               SummarizedExperiment_1.16.0
## [29] DelayedArray_0.12.0         BiocParallel_1.20.0        
## [31] matrixStats_0.55.0          Biobase_2.46.0             
## [33] GenomicRanges_1.38.0        GenomeInfoDb_1.22.0        
## [35] IRanges_2.20.1              S4Vectors_0.24.0           
## [37] BiocGenerics_0.32.0         data.table_1.12.6          
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.4-1       rprojroot_1.3-2        htmlTable_1.13.2      
##  [4] XVector_0.26.0         fs_1.3.1               base64enc_0.1-3       
##  [7] rstudioapi_0.10        farver_2.0.1           affyio_1.56.0         
## [10] bit64_0.9-7            AnnotationDbi_1.48.0   lubridate_1.7.4       
## [13] xml2_1.2.2             splines_3.6.1          geneplotter_1.64.0    
## [16] knitr_1.26             zeallot_0.1.0          Formula_1.2-3         
## [19] broom_0.5.2            annotate_1.64.0        cluster_2.1.0         
## [22] dbplyr_1.4.2           BiocManager_1.30.10    compiler_3.6.1        
## [25] httr_1.4.1             backports_1.1.5        assertthat_0.2.1      
## [28] Matrix_1.2-18          lazyeval_0.2.2         cli_1.1.0             
## [31] formatR_1.7            acepack_1.4.1          htmltools_0.4.0       
## [34] tools_3.6.1            affy_1.64.0            gtable_0.3.0          
## [37] glue_1.3.1             GenomeInfoDbData_1.2.2 Rcpp_1.0.3            
## [40] cellranger_1.1.0       vctrs_0.2.0            preprocessCore_1.48.0 
## [43] gdata_2.18.0           nlme_3.1-142           xfun_0.11             
## [46] rvest_0.3.5            testthat_2.3.0         lifecycle_0.1.0       
## [49] gtools_3.8.1           XML_3.98-1.20          zlibbioc_1.32.0       
## [52] scales_1.1.0           hms_0.5.2              lambda.r_1.2.4        
## [55] curl_4.2               yaml_2.2.0             memoise_1.1.0         
## [58] gridExtra_2.3          rpart_4.1-15           latticeExtra_0.6-28   
## [61] stringi_1.4.3          RSQLite_2.1.2          highr_0.8             
## [64] genefilter_1.68.0      checkmate_1.9.4        caTools_1.17.1.2      
## [67] rlang_0.4.2            pkgconfig_2.0.3        bitops_1.0-6          
## [70] evaluate_0.14          labeling_0.3           htmlwidgets_1.5.1     
## [73] bit_1.1-14             tidyselect_0.2.5       R6_2.4.1              
## [76] generics_0.0.2         Hmisc_4.3-0            DBI_1.0.0             
## [79] pillar_1.4.2           haven_2.2.0            foreign_0.8-72        
## [82] withr_2.1.2            survival_3.1-7         RCurl_1.95-4.12       
## [85] nnet_7.3-12            modelr_0.1.5           crayon_1.3.4          
## [88] futile.options_1.0.1   KernSmooth_2.23-16     rmarkdown_1.18        
## [91] readxl_1.3.1           locfit_1.5-9.1         blob_1.2.0            
## [94] reprex_0.3.0           digest_0.6.23          xtable_1.8-4          
## [97] munsell_0.5.0          viridisLite_0.3.0

DEG analysis, GO analysis and bubble plot preparation

Thomas Dobrenel and Nicolas Delhomme

2020-01-17

1 Setup

2 Expression data

2.1 Loading of Fold changes and adjusted p-values

2.2 Apply filters for the DEGs

2.2.1 For the 6hrs timepoint

2.2.2 For the 24hrs timepoint

2.2.3 For Wouter’s data

3 GO analysis

3.1 For the 6hr timepoint

3.2 For the 24hr timepoint

3.3 For Wouter’s data

3.4 Add the experiment name in the GO tables

4 Combine the GO tables

4.1 Prepare the tables

4.2 Keep only the biological processes

4.3 Combining Up and Down-regulated genes together

4.4 Removing the root terms

5 Create a list of interesting GO

5.1 List of GO groups overrepresented among genes induced after 24hrs of Suc starvation

5.2 List of GO groups overrepresented under the different conditions

6 Create a sublist of GO with only the GOs of interest

7 Make the bubble plot

8 Session Info